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Abstract 


A new p-version finite element formulation for steady, incompressible fluid flow and 
convective heat transfer problems is presented. The steady-state residual equations are 
obtained by considering a limiting case of the least-squares formulation for the transient 
problem. The method circumvents the Babu§ka-Brezzi condition, permitting the use of 
equal-order interpolation for velocity and pressure, without requiring the use of arbitrary 
parameters. Numerical results are presented to demonstrate the accuracy and generality 
of the method. 



1. Introduction 


Methods used to improve the quality of finite element solutions include h-methods, 
p-methods and h-p methods [l]. Finite element programs generally employ the h- version 
of the method, where the accuracy of the solution is improved by refining the mesh while 
using fixed low order element interpolation. This approach can be contrasted with spectral 
methods, which use high order global approximation functions. Spectral methods have 
been used to obtain highly accurate solutions to a variety of fluid dynamics problems [2]. 

During the past two decades p- and h-p versions of the finite element method have 
been developed which combine the geometric flexibility of standard low-order finite element 
techniques with the rapid convergence properties of spectral methods. In the p-version of 
the method, accuracy is achieved by increasing the element interpolation rather than refin- 
ing the mesh. The h-p version involves a combination of mesh and polynomial refinement. 
These finite element developments have been primarily applied to structural mechanics 
analysis. 

In recent years, ’spectral element’ methods which are similar to the p- and h-p version 
finite element methods have been developed specifically for incompressible fluid flow[3]. 
Adaptive mesh strategies for the spectral element method have recently been investi- 
gated^]. Compressible flow problems have already been solved using h-p adaptive finite 
element methods [5]. 

Finite element methods are typically based on either variational or weighted resid- 
ual formulations. Variational formulations are generally viewed as producing the ’best’ 
approximation to the exact solution of a variational problem. Unfortunately, variational 
principles cannot be constructed for the Navier-Stokes equations. In the absence of a varia- 
tional principle, a weighted residual method is usually employed, the Galerkin formulation 
being the most popular [6]. 

A standard Galerkin formulation for the steady, incompressible Navier-Stokes equa- 
tions in primitive variable form unfortunately requires the use of velocity and pressure 
interpolation functions which satisfy the Babuska-Brezzi condition. Several methods of 
circumventing this problem have been investigated and reported in the literature[7-9]. The 
approach derived here is based on the formulation presented by DeSampaio[10], although 
the present method does not require an ’upwinding’ parameter. In fact, the numerical 
studies to be presented indicate that the use of ’upwinding’ (Petrov-Galerkin weighting) 
may actually degrade the solution when high order approximation functions are used. 

The paper is org aniz ed as follows. The hierarchical p-version approximation functions 
are described in the next section. Using these p-version approximation functions, the 
finite element method is applied to the one- dimensional Burgers’ equation in section 3 
and the two-dimensional convection- conduction equation in section 4. The steady-state 
residual equations are obtained by considering the steady-state limit of the least-squares 
formulation for the transient problem. Numerical examples with known solutions are used 
to evaluate the accuracy of the method. The method is extended to the incompressible 
Navier-Stokes equations in section 5. Accurate numerical results are obtained in section 
6 using equal-order interpolation for velocity and pressure without requiring the use of 
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arbitrary parameters. Conclusions are provided in the last section of the paper. 


2. p- Version Finite Element Approximation 

The p- version approximation functions and nodal v ari ables for the two-dimensional 
element shown in figure 1 can be obtained by taking products of one- dimensional approxi- 
mation functions and nodal variable operators in the £ and 77 natural coordinate directions. 



□ Nodes with hierarchical degrees of freedom 
O Non-hierarchical nodes 


Figure 1: p-Version Element in natural coordinate space (£,77) 


The one- dimensional approximation can be obtained from a Taylor series expansion, 
as shown in Appendix 1. In the £ direction the approximation is given by: 
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Similarly, in the 77 direction: 
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Note that the first three terms on the r.h.s. of eqns (1) and (4) are the standard 
Lagrange interpolation functions for a parabolic element. Higher order terms are included 
by simply adding hierarchical degrees of freedom at the center node. 

The two-dimensional approximation functions and nodal variables can be obtained by 
simply taking products of expressions (1) and (4), which yields the element approximation: 


©*({,!)) = [iV]{0> (7) 

The contents of [N] and {9} are given in table 1. 

Adaptive p-level refinement is extremely easy when hierarchical approximation func- 
tions are used. The element nodal configuration and geometry do not change. The p-level 
is increased (or decreased) by simply adding (or subtracting) degrees of freedom at nodes 
2, 4, 6, 8 and 9 shown in figure 1. Inter-element continuity conditions are satisfied by con- 
straining appropriate degrees of freedom at the mid- side nodes 2,4,6 and 8 . 


3. Burgers’ Equation 

To illustrate the basic ideas, let’s first consider the one- dimensional Burgers’ equation: 
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Table 1: Hierarchical approximation functions and nodal variables 


Node 

Approximation 

Hierarchical Nodal 

Order of 

Number 

Functions 

Variables 

Derivatives 

1 


(0)j = 0|^=_1,,=-1 

0 

3 

N J’JVj* 

(0) 3 = 0|f=+l,7;=-l 

0 

5 

nI” jv 3 °< 

(0) 5 = €>le=+i,7 7 =+i 

0 

7 


(0)t = ©lf=-l,i»=+l 

0 


JVf’JVj* 

(0) 2 = 0| f=o ,i7=-l 

0 

2 


(® S ‘) 2 = h^) £ „ 

i=3,4,..., p 


N° 3 'N° 2 t 

(0) 6 = 0|f=O,J7=+l 

0 

6 


( 0 fOe = s(^)*«o,„-+i 

i— 3,4,..., p 



(0) 4 = 0|f= + l,7J=O 

0 

4 

NpNi* 


j=3,4,..., p 


N}'N** 

(0) 8 - 0|(— — 1,TJ = 0 

0 

8 



j=3,4,..., p 



(0)9 = ®|^=0 

0 


n^n? 

( 0 f)9 = l(w ) (=v=0 

i=3,4,..., p 

9 


(0)^)9 = j(S£)c= n =0 

j=3,4,..., p 


N^N? 

( 0 e* v )o = A [^(^)] {=J7=0 

1=3,4,... ,p 
j=3,4,...,p 
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( 8 ) 


du du d 2 u 

dt U dx V dx 2 


= 0 


Integrating with, respect to time: 



( 9 ) 


Using the trapezoidal rule to approximate the integral in eqn (9): 
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The least- squares criterion is now applied by minimizing the integral of the squared 
residual with respect to the degrees of freedom at time level ’n+l’. The resulting equations 
can be interpreted as a weighted residual formulation with weight functions as: 


{ W } = [N] : 


At 
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where 
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Note that the product rule of differentiation has been used to obtain the above ex- 
pression. In the steady-state limit, the weighted residual formulation becomes 


The above expression is clearly a Petrov-Galerkin formulation where velocity depen- 
dent ’least- squares’ type terms have been included in the weighting functions given by eqn 
(11). The time step At serves as the upwinding parameter. 

The Galerkin weighting [N] T can be used in the standard fashion to integrate the 
diffusive term by parts: 
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— — dfi -f- v I 
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\ dN ] 
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. dx . 


T du At 

dx dn+ T 



du 

dx 


r 


(13) 


Note that the boundary term only applies to portions of the boundary where the 
primary variable ’u’ has not been prescribed. 

For a one- dimensional convection- diffusion problem with a constant convective veloc- 
ity, the parameter A t may be selected such that nodally exact steady- state solutions are 
obtained. However, for general problems in fluid dynamics and heat transfer the optimal 
value of the upwinding parameter is not known. 
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Numerical Example 


Consider the solution of the steady-state Burgers’ equation with the boundary con- 
ditions u(0)=l, u(l)=0 and let v — 0.01. Solutions to equation (13) were obtained using 
Newton-Raphson iteration. The non-symmetric matrix equation resulting from (13) was 
solved using a frontal solver[ll] modified for use with p-version elements. No difficulties 
were encountered in obtaining converged solutions. 



Figure 2: Solutions to Burgers’ equation at p=2 


A unif orm mesh of ten elements at a p-level of p=2 was used to produce the oscillatory 
results shown in figure 2. The standard Galerkin formulation corresponds to the At = 0.000 
solution, which is quite oscillatory. Petrov- Galerkin weighting is commonly prescribed as a 
cure for this oscillation problem[8], and it is clear in this case that the upwinding parameter 
A t is in fact suppressing the oscillations. Unfortunately, it is also clear that none of the 
’solutions’ shown in figure 2 are very accurate. As has been previously pointed out by 
Gresho and Lee[12], the appropriate solution to this ’wiggle’ problem is adequate mesh 
refinement. Here the use of p-refinement rather than h-refinement will be explored. 

The exact solution to this problem[13] is given by: 


f 1 — exp[u(x — l)/i^] 1 
11 + exp[u(x — \)/v] j 


( 14 ) 


where u satisfies: 
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=— 7 = exp(— u/u) 
u + 1 

Using eqn (14) the L 2 norm: 

-u k y<iny (15) 

can be used to measure the solution error. 

Using (15), the error was computed from the solutions obtained using selective p- 
refinement. The p-level was changed in the last element (0.9 < x < 1.0) where the solution 
is changing rapidly and the p-level was fixed at p=2 for all the other elements. The solution 
error was computed using 12 point Gauss quadrature to evaluate the integral in eqn (15). 
The computed error values are plotted as a function of the At parameter in figure 3 for 
several p-levels. The results at p=2 confirm that the use of Petrov-Galerkin weighting 
can in fact improve the solution accuracy, provided the correct value of the upwinding 
parameter is selected. At p=3, the best possible At parameter improves the accuracy only 
slightly, and at p=4 the best parameter choice is At = 0. 

0.08 
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§ 0.04 
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0.02 


0.00 

0.00 0.02 0.04 0.06 0.08 0.10 

delta t 

Figure 3: Burgers’ equation solution error 

For this problem it appears that the standard Galerkin formulation is best when higher 
order approximations are used. In figure 4, Galerkin solutions in the oscillatory region are 
compared with the exact solution. The oscillations are quite small at p=4, and at p=8 the 
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solution is essentially exact. Highly accurate solutions have been obtained without the use 
of upwinding by simply employing local p-refinement. 



Figure 4: Solutions to Burgers’ equation at higher p-levels 


4. Convection-Conduction Equation 


Before proceeding to the Navier-Stokes equations, let’s first examine the simpler 
convection-conduction equation in two dimensions: 


dT dT dT 1 (d 2 T d 2 T\ 
dt U dx V dy Pe\ dx 2 dy 2 ) 


(16) 


where u and v are known dimensionless velocity components and T is the dimensionless 
temperature. Applying the same approach to eqn (16) as was applied to Burgers’ equation 
(8) yields the following steady-state result: 


g x J ax l dy J dy) 


/ 0 M'(-I + ’S)" + s:/ 0 ([ 

+ f L{*} TdSl = TeJ r [ N } TvT - ndr 

here n is the outward normal vector at the boundary, r is the residual: 


(17) 
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3T h dT h 1 (d 2 T h d 2 T h \ 
T ~~ U ~dx + V ~dy ~ Pe V dx 2 + dy 2 ) 

and the vector {$} is given by: 


f , 1 dr 

\dN l 

T 

dN- \T 

1 

rd 2 Ni T 

1 

d 2 Ni 

l $ j “ d{T} ~ u 

. dx . 

+ v 

■ dy J 

_ Te 

. dx 2 J 

~Te 

dy 2 J 


Once again, the boundary term in (17) only needs to be computed on portions of the 
boundary where the primary variable ’T’ is not specified. 

Numerical Example 

The well known Smith-Hutton problem[14] will be solved using eight uniform p-version 
elements, shown in figure 5. The condition at the inlet is prescribed as: 

T = 1 + tanh[10(2x + 1)]; y = 0, -1 < x < 0 (inlet) (18) 



The dependent variable ’T’ is essentially zero on x = ±1 and y = 1 and very nearly 2 
at the origin of the coordinates. The climb from 0 to 2 occurs very sharply halfway along 
the inlet. A zero flux boundary condition (q = VT ■ n = 0)is assumed at the outlet. For 
this problem the velocity field is specified analytically as: 

u = 2y(l — x 2 ), v = — 2x(l — y 2 ) (19) 
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corresponding to a streamfunction 


# = -(l-z 2 )(l-j, 2 ) 


( 20 ) 


In the case of pure convection (Pe — * oo), the exact solution to this problem is easily 
obtained in terms of the streamfunction 'J': 


T = T(¥) = 1 + <anh[10(l - 2\/l + ¥)] 


( 21 ) 


For very large Peclet numbers 
and the L 2 norm 


(Pe > 10 6 ), eqn (21) can be used as an exact solution 


Hello 



( 22 ) 


can be used to measure the solution error. 

It is well known that the standard Galerkin formulation is optimal for pure conduction 
(Pe — > 0) problems. It is commonly held that the Galerkin formulation is acceptable in 
low Pe cases, but inadequate for convection dominated problems. The availability of an 
exact solution for the convection dominated case permits an objective evaluation of the 
optimal upwinding parameter to be used with p- version finite elements. 



Figure 6: Smith-Hutton problem error 
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delta t = 0.00 
delta t = 0.01 
delta t = 0.10 


Figure 7: Solutions at Pe 



- FEM solutions 

• Pe=1 0 reference 
■ Pe=1 00 reference 

♦ Pe=500 reference 
▼ Pe=1 000 reference 
a Pe=1 0 A 6 reference 


Figure 8: Smith-Hutton problem solutions 


Eqn (17) was placed in matrix form and solved using the same p- version non- symmetric 
matrix solver used in section 3. Solutions obtained with Pe = 10 6 were used in eqn (22) 





to compute the error for several values of the upwinding parameter At. The results for 
several p-levels are presented in figure 6. Results for p < 4 are not presented because the 
inlet boundary conditions cannot be properly represented at low p-levels with the present 
mesh. From figure 6 it is clear that the use of Petrov-Galerkin weighting increases the 
solution error. The standard Galerkin formulation (At = 0.0) produces the best results at 
these p-levels. 

The effect of the use of Petrov-Galerkin weighting at p=6 is illustrated in figure 7, 
where the upwinding parameter At is seen to introduce rather than eliminate oscillations 
in the outlet solution. Petrov-Galerkin weighting was also tested at a low Peclet number 
(Pe=10) and found to be of no benefit. 

Galerkin solutions at the outlet for a range of Peclet numbers are presented in figure 8 
along with reference solutions [14] considered to be reasonably accurate. The finite element 
solutions presented were obtained with p=6 for all cases except Pe=1000, where it was 
necessary to use p=7. It should be noted that these results are quite superior to previously 
published p- version least-squares results[15j. 


5. Incompressible Fluid Flow 


The p- version finite element method used in sections 3 and 4 will now be extended to 
the incompressible Navier-Stokes equations. Two-dimensional, incompressible fluid flow 
with constant properties can be described by the following set of equations expressed in 
dimensionless form: 


du du du dP 1 fd 2 u 9 2 u\ _ 

dt + U dx dy ^ dx Re\dx 2 + dy 2 ) 

dv dv dv dP 1 / d 2 v _ 

dt U dx + V dy dy Re V dx 2 ^ dy 2 ) 


du 

dx 


dv 

+ = 0 

dy 


(23 a) 
(236) 
(24) 


The approach used in sections 3 and 4 can be applied to the momentum equations 
(23) to give: 


u n+1 ~u n +Y ( f * + / i n+1 ) = 0 (25a) 

v n+i -v n +^ (/ 2 n + / 2 n+1 ) = 0 (256) 

where 


h 


du du dP 1 /d 2 u d 2 u\ 

U dx + V dy + dx Re V dx 2 dy 2 ) 


( 26 a) 
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(266) 


dv dv dP 1 fd 2 v d 2 v\ 

2 dx dy + dy Re V dx 2 dy 2 / 

Let 7 ’i, 7 , 2 and r 3 represent the residuals in equations (25a), (25b) and (24), respectively. 
Applying the least squares criterion at time level ’n-f-1’ to the integrated sum of the squared 
residuals yields: 


Li 


dry dr 2 dr ^ 

T 1 + a r ~r r 2 + -zrr^ r 3 


n ld{u}' J ' d{u} ' d{u} 
dr i dr 2 dr. 


<m = {0} 


f dry dr 2 dr z 

/ ar"T ri + aTT r2 + 57~\ r 3 dU = {0} 

J u id{v} d{v } d{v} J 


In 


F dry dr 2 dr z 

^*1 • or m ^2 ( ^*3 


]<ffi = {0} 


(27a) 

(276) 

(27c) 


.d{py i ' 5{p> " d{p} 

Expressions for the derivatives of the residuals are obtained assuming equal order 
interpolation for field variables u,v and P. 

From eqn (25a): 


^ =[jvr + Aia/ “ +I 


d{u } l ”' ■ 2 a{ti} ’ 

Similarly, from eqn (25b): 
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And from eqn (24): 

dr 3 


= [NY + 


d{u} 


r dN 
dx J 


dr 3 


r3 _ r 5 ^i T _ f n \ 

“ l dv J 5 dm “ l l 


d{v} L dy 


dr 2 


At d/ 2 n+1 


2 d{v} ’ d{P} 2 d{P} 


dr 3 
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(286) 


(28c) 


Substituting eqns (28) into (27) and considering the steady-state limit yields the 
following system of equations: 


f / r , T At df\ \ At d/ 2 , \dNT, T rdu dv\ r 1 . . 

/JO') ^Tsra)'-*T3mMd (s*s)r”"W <”•) 


B f>'- * K*tS» * [f ]'(** 5)]'“- W « 
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( 30 ) 



du du 

u-x- + v— + 
dx oy 


dP 

dx 


dv dv dP 

u ai +v di + ~di 


1 / d 2 u d 2 u 
Re vdx 2 + dy 2 


dPl -f- 


1 /d*v 
Re V dx 2 + 



where the factor ^ has been cancelled out of eqn (30). 

The above derivation is similar to the formulation presented by DeSampaio[10]. How- 
ever, DeSampaio did not include the continuity equation (24) directly in the least squares 
minimization. As a result, DeSampaio’s formulation requires the selection of a nonzero 
value of At in order to circumvent the Babuska-Brezzi condition. Because the At param- 
eter was not helpful in the p- version finite element solutions of Burgers’ equation (section 
3) and the convection- conduction equation (section 4) it is unlikely to be of benefit here. 
With At set to zero, equations (29) represent a Galerkin treatment of the momentum equa- 
tions modified by a least-squares term to enforce the continuity equation. Equation (30) 
represents a least-squares minimization with respect to the pressure degrees of freedom. 
Setting At to zero and integrating the viscous terms in equations (29) by parts yields: 


f r 1 T / du du dP\ jr . 1 t/rdN^du \dN-\ T du\ 

/ n H v‘di + '’ai + 9^) dil + Tcj a \ bd ai + la^J ai> 


f\dNi T fdu dv\^ 1 f \„1 T „ ~ 

+ /Jld yVx + dy) dQ= TeJ r H Vu ‘ n<fr 

/ r * t 1 r / dv dv dP\ 1 f (\ dNi T dv rdlV]- 1 Ov\ 

J n [ N \ { u di +v fr, + fa) dSl + TeL { lad ai + bd ad 


+ 


9x dy dy ) Re V l dx J dx L dy J dy 

f\ dN^/du dv\ jr ^ 1 f r JT , 

/Jad (ai + ay) d{l = Rij r r J Vv ' ndr 


(31a) 


(316) 


where n is the outward normal vector at the boundary. Equations 30, 31a and 31b consti- 
tute the required set of equations. 


6. Numerical Examples 

Two numerical examples are presented in this section to demonstrate the accuracy 
and convergence characteristics of this p- version formulation. 

a) Example 1: Driven Cavity, Re=1000 

Consider the well known ’driven cavity’ problem shown in figure 9. The velocity 
components are constrained along the entire boundary and the pressure is constrained to 
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The coupled, nonlinear system of equations (30, 31a and 31b) was solved by succcessive 
substitution using the p- version finite element non-symmetric matrix frontal solver used 
in the previous sections. The p-levels were sequentially increased starting with p=2. A 
starting vector of {£} = {0} was assumed at p = 2. For p > 3, the solution obtained 
at the previous p-level was used to generate an initial solution vector. This procedure 
is quite easy to implement because of the hierarchical nature of the degrees of freedom. 
The solution was considered to be sufficiently converged (in the iterative sense) when the 
maximum change in the solution vector was less than 10 -4 . Generally around 10 iterations 
were required at each p-level. 



X 

Figure 11: Plots of v(x,y=0.5) for driven cavity 


The solutions for ’u’ along the vertical centerline (x=0.5) are presented in figure 10 
and the solutions for ’v’ along the horizontal centerline are shown in figure 11. The results 
reported by Ghia et.al.[16] are also presented as a reference. Note that the solutions 
converge rapidly and the agreement with the reference solution is excellent. 

The vector and streamline plot is shown in figure 12 and the pressure contours are pre- 
sented in figure 13. Both the velocity and pressure solutions are quite smooth, even though 
equal order interpolation has been used without upwinding or other special techniques. 

b) Example 2: Asymmetric Sudden Expansion, Re=229 

In the second example a 3:2 asymmetric sudden expansion problem is examined 
and the numerical results are compared with the experimental results of Denham and 
Patrick[17]. The highest Reynolds number from the experiment, Re=229, is considered, 
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Figure 13: Pressure contour plot for driven cavity 
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where the Reynolds number is based on the mean velocity and channel half-width at 
the inlet. This problem was selected because the experimental data is very nearly two- 
dimensional and steady in the expansion region. At higher Reynolds numbers, the flow 
pattern is generally three-dimensional and often not truly laminar. 

The problem geometry (expressed in step height units), boundary conditions and finite 
element mesh are shown in figure 14. The inlet boundary conditions for ’u’ were determined 
using a least-squares fit to the inlet velocity measurements, and the vertical velocity V 
was set to zero. Fully developed conditions were imposed 40 step heights downstream. 
In addition to the natural boundary conditions (Vu. • n — V v • n — 0), the pressure was 
constrained to a reference value of P — 0 at the outlet. 


u=0 v=0 

P=0 

8u/3x=0 
dv!dx=0 

5 10 U=° v=0 20 40 

Figure 14: Asymmetric Sudden Expansion Problem 


u=f(y) 

v=0 



Figure 15: Plot of viscous shear stress along lower wall 


Computations were performed up to a p-level of p=7 using the same iterative pro- 
cedure described in the first example. In figure 15 the dimensionless viscous shear stress 
[r xy = + fj) along the lower wall (y=-l) is presented for several p-levels. The converged 
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Figure 16: Velocity Profiles 


solution indicates a recirculation zone length L r prediction of L r = 9.3, which agrees fairly 
well with the experimental value ( L r ~ 9.8). 

The converged velocity profiles at x=0,2,4,6 and 8 are compared with the experimental 
results in figure 16. The agreement with the experimental values is excellent for x=0 
through x=6 and good for the x=8 profile. Similar numerical results were obtained using 
a p-version least-squares formulation [18]. Some error (~ 2%) may have been introduced 
in the process of extracting the velocity values from the graphical results of Denham and 
Patrick. 



Figure 17: Streamlines for the sudden expansion problem 


The recirculation zone streamline plot is shown in figure 17 and the pressure contours 
are presented in figure 18. Again, smooth solutions have been obtained using equal order 
interpolation without requiring the use of arbitrary parameters. 
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Figure 18: Pressure contours for the sudden expansion problem 


7 . Conclusions 

Petrov- Galer kin formulations for the steady solution of Burgers’ equation, convec- 
tive heat transfer, and incompressible fluid dynamics problems have been presented. The 
weighted residual formulations were obtained by applying the least- squares approach to 
the transient problems, with the temporal derivatives approximated by finite differences. 
In the steady-state limit, this approach leads to a Petrov- Galerkin statement with the time 
step At serving as the upwinding parameter. 

Numerical results using p-version finite elements for the one-dimensional Burgers’ 
equation and a two-dimensional convection- conduction problem revealed that the standard 
Galerkin method (At = 0) performed better than the Petrov- Galerkin formulation when 
higher order approximations ( p > 4) were used. This evaluation was made by comparing 
the Z <2 norm of the solution errors. 

Based on these numerical findings, a value of At — 0 was selected for the incompress- 
ible Navier-Stokes formulation, resulting in a method of circumventing the Babuska-Brezzi 
condition without requiring the selection of parameter values. The numerical solutions 
obtained using this formulation agree well with reference solutions and experimental data. 

The following specific conclusions can be made concerning these results: 

1. The Babuska-Brezzi condition can be circumvented in a very simple manner without 
the use of arbitrary parameters. 

2. ’Upwinding’(Petrov-Galerkin weighting) is only helpful if the element approxima- 
tion is not capable of accurately representing the solution. 

3. The p-version of the finite element method is well-suited to problems involving 
incompressible fluid flow and convective heat transfer, producing non-oscillatory and highly 
accurate solutions. 

4. Although h-p methods may be necessary for problems with extremely sharp gra- 
dients, the p-version of the finite element method can be very effective at increasing the 
local order of approximation to produce accurate solutions to a range of incompressible 
flow problems. 

The primary goal of this work was to demonstrate the accuracy and generality of 
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this new method. The computational efficiency of the method can be improved by using 
the technique of ’static condensation’, where the degrees of freedom on element interiors 
axe computed after determining the solution on the element boundaries. Also, for large 
scale problems, indirect methods are more efficient than direct methods such as the frontal 
solver. Future work will concentrate on these efficiency improvements and the development 
of reliable error indicators for use with adaptive methods. 
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Appendix A 


Consider a three node, one dimensional element: node 1 at £ — 1, node2 at £ 0 

and node3 at £ = +1, where £ is a local coordinate. A dependent variable 0(0 can be 
expanded in a Taylor series about node 2 (£ = 0): 
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Subtracting eqn (A. 2) from (A. 3) and rearranging gives 
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Adding eqn (A. 2) to (A. 3) and rearranging gives 
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Substituting eqns (A. 4) and (A. 5) into (A.l) yields 
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which can be rewritten as 
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where 

^ f 2 for even i 
t 1 for odd i 

If terms of up to order ’p’ are retained, then 0(£) can be approximated as: 
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